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ABSTRACT 

Magnetic field strengths inferred for relativistic outflows including gamma-ray bursts (GRB) and active 
galactic nuclei (AGN) are larger than naively expected by orders of magnitude. We present three-dimensional 
relativistic magnetohydrodynamics (MHD) simulations demonstrating amplification and saturation of magnetic 
field by a macroscopic turbulent dynamo triggered by the Kelvin-Helmholtz shear instability. We find rapid 
growth of electromagnetic energy due to the stretching and folding of field lines in the turbulent velocity field 
resulting from non-linear development of the instability. Using conditions relevant for GRB internal shocks and 
late phases of GRB afterglow, we obtain amplification of the electromagnetic energy fraction to ~ 5 x lO""*. 
This value decays slowly after the shear is dissipated and appears to be largely independent of the initial field 
strength. The conditions required for operation of the dynamo are the presence of velocity shear and some seed 
magnetization both of which are expected to be commonplace. We also find that the turbulent kinetic energy 
spectrum for the case studied obeys Kolmogorov's 5/3 law and that the electromagnetic energy spectrum is 
essentially flat with the bulk of the electromagnetic energy at small scales. 

Subject headings: instabilities - magnetic fields - MHD - methods: numerical - relativity - turbulence - 
gamma rays: bursts 



1. INTRODUCTION 

Strong magnetic fields are required in GRB and AGN out- 
flows to enable sufficient production of non-thermal radia- 
tion. Synchrotron models for GRB afterglow emission require 
magnetic field strengths many orders of magnitude larger 
than expected in smoothly shock compressed circum-burst 
medium. For the prompt GRB emission, magnetic field may 
be advected from the central engine to the emission region, 
but the degree of advected magnetization is uncertain and 
may be small. The source of sufficiently strong magnetic 
field responsible for non-thermal GRB emission has thus re- 
mained a mystery, and a similar situation exists for AGN. 
Possible mechanisms capable of generating magnetic fields in 
relativistic collisionless shocks are plasma instabilities, such 
as th e Weibel two-stream inst ability (Gruzinov & Waxman 
119991: iMedvedev & Loeblll999l) . Rece nt plasma simula tions 
using the particle-in-cell (PIC) method (lBunemanlll993h have 
demonstrated ma gnetic field generation in relativistic col - 
Hsionless shocks (INishikawa et al] 120031: ISpitkovskvll2008l) . 
However, the size of the simulation boxes is orders of mag- 
nitude smaller than the GRB and AGN emitting regions. It 
remains unclear whether fields generated on scales of tens of 
plasma skin depths, as in the current PIC simulations, will 
persist at sufficient strength, or at all, on radiation emission 
scales. In addition, the recent plasma simulations ( Spitkovskv. 
|2008|) are two-dimensional hence missing intrinsically three- 
dimensi onal effects whi ch may be necessary for dynamo op- 
eration (lGruzinovll2008h . Thus it is far from clear that mag- 
netized regions with size scales and magnitudes inferred for 
AGN and GRBs can be generated by plasma instabilities 
alone. 

Another possibility is that weak magnetic field is 
greatly amplified by a magnetohydrodynamic (MHD) tur- 
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bulent dynamo (j Kazantsev 1968: Kulsrud & Andersonlll992t 
iBoldvrev & Cattaneo.2004; .Gruzinov,2008) . In our scenario 
the turbulence is driven by non-linear development of the 
Kelvin-Helmholtz instability present due to velocity shear 
Shear is naturally expected in realistic explosive outflows due 
to intermittency and inhomogeneity in the outflow itself and in 
the surrounding medium. Collisions in GRB internal shocks, 
for example, likely occur between blobs of differing size and 
shape resulting in oblique shocks and vorticity generation in 
the post-shock flow (Goodman & MacFadven 2008). In jets, 
shear is intrinsic between a fast jet core and the surrounding 
"cocoon" which is often turbulent due to Kelvin-Helmholtz 
instability. Shocks launched into a clumpy medium, as ex- 
pected, for example, in the regions surrounding massive stars, 
will lead to significant vort icity generation a nd velocity shear 
during clump shocking (Goodman & MacFadven 2008) or 
from ion streaming ( Couch, Milosavljevic & Nakar 2008). In 
these cases the Kelvin-Helmholtz instabiUty is Ukely to gen- 
erate turbulence from the shear flow. These flows are highly 
ionized and as such are expected to contain some magnetic 
field however weak. 

In this letter, we present high-resolution, three-dimensional 
numerical simulations of relativistic MHD demonstrating am- 
plification of weak magnetic fields by Kelvin-Helmholtz in- 
stability induced turbulence. This process is demonstrated to 
be viable as the origin of strong magnetic fields present behind 
relativistic collisionless shocks of astrophysical jets in AGN 
and GRBs. Our simulations also allow us to study the basic 
properties of turbulence in the largely unexplored relativistic 
MHD case. In §|2]we describe the initial setup of our simu- 
lations. In §[3] we present the simulation results. We further 
discuss the results in § H] 



2. MHD EQUATIONS AND INITIAL SETUP 

The system of special relativistic MHD (SRMHD) is 
governed by the equations of mass conservation, energy- 
momentum conservation and magnetic induction, respec- 
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Fig. 1 . — Plasma /3, the ratio of gas pressure to the magnetic pressure, for Model A- 1024. Results at ? = 2 and / = 12 are shown in the left and right panel, 
respectively. 1024^^ numerical cells are used in this calculation for which the numerical box size and light crossing time is 1. At ? = 0, the speed of the gas is set 
to Vv = 0.5 in the top and bottom quarters of the box and v,- = —0.5 in the middle half of the box, with small perturbations added. Initially, the box is filled with a 
weak magnetic field, which corresponds to /3 = 2 X 10*. By / = 2, vortices due to the Kelvin-Helmholtz instability are clearly visible. The instability grows and 
the magnetic field is amplified by stretching and folding of field lines. By ( = 12, the gas has become fully turbulent due to the instability and the plasma beta is 
highly variable with values as low as ~ 10 in some regions. 



lively given by 

V,pm" = 0, (1) 
V,.r^'' = 0, (2) 

-B = V X (v-x B). (3) 
of 

Here p is the mass density in the local rest frame, is the 
four-velocity, T^" is the energy-momentum tensor, B is the 
magnetic field three-vector in the lab frame, and v is the three- 
velocity. The energy-momentum tensor, which includes both 
the fluid and electromagnetic parts is given by, 

T^"" = phu^'u" + Pg'"' + (u'^u" + ^g'"')b^ - b^f , (4) 

where P is the pressure, g'^'^ is the metric, b^^ is the four- vector 
of the magnetic field measured in the local rest frame, and 
h=l+e + P/pis the specific enthalpy here e is the specific 
internal energy. Units in which the speed of light is set to 
c = 1 are used, and the magnetic field is redefined to absorb a 
factor of 1 / \/47r. The system is closed by an equation of state. 
Furthermore, there is the solenoidal constraint V • B = 0. 

A mildly relativistic and weakly magnetized shear flow 
with zero net momentum was set up as the initial condition. 
A three-dimensional grid {x,y,z) in Cartesian coordinates was 
used with periodic boundary conditions applied in all direc- 
tions. The dimensionless size of the box was 1 in each di- 
rection and the center of the box was at (0,0,0). Initially the 
numerical box was filled with relativistically hot gas of con- 
stant pressure of Pq = 1 and constant density of po = 1 . The 
equation of state was assumed to be an ideal gas with gamma- 
law, f = (7- l)pe, with a constant adiabatic index 7 = 4/3. A 
shearing flow initial condition was imposed with the velocity 
in the top quarter (y > 0.25) and bottom quarter (y < -0.25) 
of the box in the positive x-direction at a speed of = 0.5 and 
in the negative x-direction at a speed of = -0.5 in the mid- 
dle half (-0.25 < y < 0.25) of the box. This corresponds to 



a shear flow with relative Lorentz factor of 2.29. The system 
had an initially weak magnetic field of B^ = 10"^ and 10"^ for 
Models A and E, respectively. Small velocity perturbations 
were added of the form, 

v,i =flpCos(k-r + (/),), (5) 

where / denotes x-, y- and z-direction, r is the position and 
(pi is a random number in the range of [0,27r). The ampli- 
tude of the perturbation is set to Op = 0.01. The wave vector 
k is given by k = 27r(cos6'pef + sin^pev), here 6p is set to 0.05. 
The initial perturbation is roughly a single mode. The rel- 
ativistic sound speed in the initial models is = 0.516 and 
the flow speed is 0.5 (with perturbations of < 0.01). Thus 
the initial shear flow is transonic. Also the initial conditions 
of the shear flow correspond to (3 = 1.2 x lO*" and 1.2 x 10"* 
for Models A and E, respectively. Here the plasma beta is 
/3 = P/Pm, where P is gas pressure and P,„ is magnetic pres- 
sure. In the absence of magnetic field, Kelvin-Helmholtz in- 
stability develops no matter how small or large the velocity 
shear is, whereas strong magnetic field can suppress the insta- 
bility. This is true in 3D even for supersonic flow speeds be- 
cause there are always d irections al ong which the shear flow 
component is subsonic jLandau &^ Lifshitz 1959). Because 
the initial field in our models is small and dynamically unim- 
portant, Kelvin-Helmholtz instability develops in our 3D sim- 
ulations. 

3. RESULTS 

The SRMHD simulations in this letter were performed with 
the general relativistic MHD code Nvwa (Zhang & Wang in 
prep). The three-dimensional numerical simulations were run 
with various resolutions: 512^, 768^* and 1024-' numerical 
cells. We use A-n and E-n to denote Models A and E, re- 
spectively, where n stands for the number of numerical cells 
in each dimension. For example. Model A- 1024 stands for 
Model A with a resolution of 1024^ numerical cefls. 
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Fig. 2. — Ratio of electromagnetic and kinetic energy to total energy as 
a function of time for Models A and E. The lines for kinetic energy start at 
e ~ 4 X 10"^ at / = for all models, whereas the lines for electromagnetic 
energy start at e ~ 10"' and 10"' for Models A and E respectively. Different 
lines are for different resolutions: 1024^^ cells (solid line), 768^ cells (dashed 
lines), and 512^ cells (dotted lines). Models A and E are shown in thick and 
thin Hues, respectively. The ratio of electromagnetic energy to total energy 
(ta) grows due to growth of the Kelvin-Helmholtz instability and ensuing 
turbulence. The electromagnetic energy saturates at f ~ 9 depending upon the 
model and resolution. The kinetic energy decays over time, and it becomes 
comparable to electromagnetic energy at / ~ 14 for Model A-1024. 

3.1. Amplification of Magnetic Fields 

In all models the initial discontinuity in velocity triggered 
Kelvin-Helmholtz instability. Two snapshots of Model A- 
1024 are shown in Fig. [T] Kelvin-Helmholtz vortices start 
to appear on the plane of the velocity discontinuity after ~ 2 
light crossing times of the numerical box due to the roll-up of 
the shear flow. The magnetic field at the shear layer is being 
folded and twisted by the vortical fluid motion. The strength 
of the magnetic field starts to increase rapidly, growing ex- 
ponentially as shown in Fig. [2] As the instability evolves, 
more and more gas is turned into a state of turbulence. At 
about t ^ 7, the whole box has become fully turbulent, and 
the electromagnetic energy starts to saturate. The turbulence 
eventually decays because there is no bulk shear flow left to 
further drive the turbulence. Fig.|2]shows the evolution of the 
ratio of the electromagnetic energy to total energy for Mod- 
els A and E at various levels of resolution, here the total en- 
ergy does not include rest mass energy. Also shown in Fig.|2] 
is the ratio of the kinetic energy^ to total energy. Note that 
the kinetic energy decays faster than the electromagnetic en- 
ergy. For Model A- 1 024, the electromagnetic energy becomes 
comparable to kinetic energy at f 14. 

The spatial structure of plasma beta for Model A-1024 is 
shown in Fig[T]at two times. The structure is clumpy and fila- 
mentary as expected for magnetic fields amplified by a turbu- 
lent dynamo. At the stage of full turbulence, the average ratio 
of the electromagnetic energy to total energy is ^ 2 x lO"-' 
for Model A-1024, whereas maximum values of es.m-dx ^ 0.2 
exist throughout the box. Thus magnetic field is very strong 
in local patches of the turbulence. These strongly magne- 
tized clumps are preferably elongated along the magnetic field 
direction due to the intrinsi c anisotropy of MHD turbulence 
dGoldreich & Sridhai|[l995h . It should be noted that there is 
no overall magnetic flux except for a very low magnetic flux 

^ The kinetic energy density is defined as pr(r— 1), where p is mass den- 
sity in the fluid rest frame and F is the Lorentz factor in the laboratory frame. 
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Fig. 3. — Spherically-integrated spectrum of (a) kinetic energy (top) and (b) 
electromagnetic energy (middle), and (c) ratio of electromagnetic spectrum 
to kinetic energy spectrum (bottom) for Model A-1024. Different lines are 
for different times: t = 4 (solid lines), f = 6 (long dashed lines), f = 8 (short 
dashed lines), f = 10 (dotted lines), t = 12 (short dash-dotted lines), and ^ = 14 
(long dash-dotted lines). Also shown in (a) is a line representing the power 
law E(k) ~ IC^f^ for comparison. 

in jc-direction present in the initial conditions. The strength of 
the magnetic field is a result of folding and twisting of existing 
field lines by turbulent motions. This process increases mag- 
netic energy density without increasing total magnetic flux. 

To study the effects of numerical resolution and initial mag- 
netic field strength on the level of magnetization that can be 
reached, we have run a series of calculations. Models A-768, 
A-512, E-768 and E-512 for comparison (Fig. |2]i. As ex- 
pected, the onset of the instability is delayed and the satu- 
rated level of magnetization is lower for lower resolution runs 
due to their higher numerical diffusivity. The maximal over- 
all ratio of the electromagnetic energy to total energy reaches 
8.2 X lO"-*, 1.6 X 10"^^ and 2.4 x 10"^^ for Models A-512, A- 
768 and A-1024, respectively. For Model E, which started 
with higher initial magnetic field, the final magnetization is 
higher than that of Model A as we also expected because less 
stretching and folding is required for reaching the same mag- 
netic field strength. The maximal overall ratio of the elec- 
tromagnetic energy to total energy reaches 4.1 x lO""* and 
4.5 X 10"^ for Models E-512 and E-768, respectively. But 
it should be noted that the two orders of magnitude difference 
in the initial fields causes only a factor of < 3 in the final ratio 
of the electromagnetic energy to total energy. These results 
suggest that the "true" answer is ~ 5 x lO"-' for the chosen 
conditions. 

3.2. MHD Turbulence 

We now analyze the properties of decaying mildly relativis- 
tic adiabatic MHD turbulence using the data of Model A- 
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1024. The initial flow is mildly relativistic and transonic 
with shear velocity of 0.5 corresponding to a Lorentz factor 
of To = 1.15. During the evolution, the velocities of plasma 
elements decrease due to dissipation into heat and the work 
done to stretch and fold magnetic field lines. The typical ve- 
locities in Model A-1024 at f = 4, 8, 14, are about 0.4, 0.3 
and 0.1, with a standard deviation of 0.1, 0.1 and 0.05 re- 
spectively. Hence, the MHD turbulence in our simulations is 
mostly subsonic except for at the very beginning. 
We define the kinetic energy spectrum Ek(k) as 

j Ek{k)dk={pT(T-\)), (6) 

where F is the Lorentz factor, and the electromagnetic energy 
spectrum E^ik) as 

j E„{k)dk={TZ). (7) 

where is the 00-componentof the energy-momentum ten- 
sor for electromagnetic field. 

It is clearly shown in Fig.[3]that the kinetic energy spectra 
follow the Kolmogorov's 5/3 law in the inertial range and fall 
off at k/li: ^ 100. The spectrum is slightly flatter than k~^l^ 
before the saturation of the overall electromagnetic magnetic 
energy at f ~ 9, whereas it is slightly steeper than k~^l^ after 
the saturation. The kinetic energy decreases over time due to 
viscous dissipation. 

The electromagnetic energy spectrum is also shown in 
Fig. 13 Note that the electromagnetic energy spectrum in- 
creases over time until it saturates. But the overall shapes 
of the spectrum at f = 4 and 6 are very similar. Th is is a typ- 
ical b ehavior of the so-called small-scale dynamo (iKazantsevI 
Il968l) . The electromagnetic energy saturates at small scales 
ik/2-K > 100) and large scales {k/liT < 5) first, then at in- 
termediate scales (5 < k/ln < 100) later. There appears to 
be both a forward cascade of electromagnetic energy from 
large scales to intermediate scales and an inverse cascade from 
small scales to intermediate scales. It is striking that the elec- 
tromagnetic energy spectra are flat and evolve very slowly af- 
ter f = 8, the time when the electromagnetic energy saturates. 
It is not surprising that f = 8 is also the time when the whole 
simulation box is filled with turbulence. A flat magnetic en- 
ergy spect rum is commonly seen in turbulent dynamo si mula- 
tions (e.g.. lBrandenburgl2001l:ISchekochihin et al.l2004 . The 
shape of electromagnetic spectra clearly indicate that the bulk 
of the electromagnetic energy is at small scales. 

Fig.|3]shows that the ratio of electromagnetic to kinetic en- 
ergy spectrum increases monotonically on almost all scales as 
the instability and turbulence develop. The electromagnetic 
energy dominates kinetic energy at small scales, whereas the 
kinetic energy dominates at large scales. The equipartition 
point moves toward larger scales during the evolution. Both 
electromagnetic and kinetic energy are decreasing after f ~ 9 
(Fig. |2|i. However, the kinetic energy decreases much faster 
than the electromagnetic energy. This means that the electro- 
magnetic resistive dissipation is less efficient than the viscous 
dissipation of kinetic energy. Therefore the magnetic fields 
amplified by the turbulent dynamo can persist for a very long 
time. 

4. DISCUSSION 



In this letter, we have presented a series of relativistic MHD 
simulations with an initial density of po = 1, pressure of fo = 1^ 
and shear velocity of vo = 0.5. The simulations are represen- 
tative of conditions in which a turbulent dynamo operates be- 
hind a relativistic astrophysical shock. The initial conditions 
for models A and E presented here are relevant for GRB inter- 
nal shocks, the late stage of GRB afterglows, trans-relativistic 
explosi ons of 1998-b w like supernovae, aka "hypernovae" 
( Kulka mTet al.lll998|). and shock breakou t from Type Ibc su- 
pernova SN2008 ("Soderberg et al."2008'). Our calculations 
predict the value of the ratio of magnetic energy to total en- 
ergy, eg = 5 X 10"^, which is typically treated as a free pa- 
rameter in radiation modelling. Further work is under way to 
cover more parameter space relevant to additional stages of 
GRB and AGN outflows, and other astrophysical relativistic 
outflows such as t hose from magnetars (Gelfand et al. 200^ 
and microquasars jMirabel et al.lll99^ . These results and a 
comprehensive study of relativistic turbulence now underway 
will be presented in future publications. 

Our simulations show that the electromagnetic energy 
structure is very clumpy. This turbulent shearing gas with 
moving magnetic clouds is likely a site for Fermi acceleration 
of charged particles, and thus a source of high energy cosmic 
rays. We are currently calculating charged particle accelera- 
tion in model A (Zrake et al., in prep). It is possible that a 
self-consistent model of GRB afterglows, which includes the 
acceleration of nonthermal electrons and the amplification of 
magnetic fields, can be built upon the processes considered 
in this letter. The spatial variability and large fluctuations of 
the turbulent magnetic field also have implications for non- 
thermal radiation which may have observable consequences. 
Synchrotron and inverse Compton spectra from the turbulent 
fields are being calculated and will be presented in a future 
paper (van Velzen et al., in prep). 

The turbulent dynamo studied in this letter may also take 
place in other astrophysical environments. It is still a mys- 
tery how the cosmic magnet ic fields in galaxies and c lusters 
of galaxies got started (see iKulsrud & Zweibel l2008l for a 
review). One possibility is that the cosmic magnetic fields 
are generated by turbulence during hier archical str ucture for- 
mation (Pudritz&Silk 1989; Kulsrud e t~anil997l) . Another 
possibility is that the strong magnetic field in the jets of 
AGN could spread to the entire universe (iRees & SettiilT968l: 
lDalv&Loeblll990l) . The Kelvin-Helmholtz induced fiirbu- 
lent dynamo we have considered can be at work in both 
cases to provide a seed field for the mean-field dynamo 
( IKulsrud & Zweibelll2008l and references therein). 



We are greatly indebted to Andrei Gruzinov for many 
stimulating discussions. We would also like to thank Yosi 
Gelfand and Martin Pessah for useful discussions. This re- 
search utilized resources at the New York Center for Compu- 
tational Sciences at Stony Brook University/Brookhaven Na- 
tional Laboratory which is supported by the U.S. Department 
of Energy under Contract No. DE-AC02-98CH10886 and by 
the State of New York and the CCNI, supported by the New 
York State Foundation for Science, Technology and Innova- 
tion (NYSTAR). 



5 



REFERENCES 



Boldyrev, S., & Cattaneo, R 2004, Physical Review Letters, 92, 144501 
Brandenburg, A. 2001, ApJ, 550, 824 

Buneman, O. 1993, in Computer Space Plasma Physics: Simulation 
Techniques and Software, ed. H. Matsumoto & Y. Omura (Tokyo: Terra 
Scientific Publ. Co.), 67 

Couch, S. M., Milosavljevic, M., & Nakar, E. 2008, ApJ, 688, 462 

Daly, R. A., & Loeb, A. 1990, ApJ, 364, 451 

Gelfand, J. D., et al. 2005, ApJ, 634, L89 

Goldreich, P, & Sridhar, S. 1995, ApJ, 438, 763 

Goodman, J., & MacFadyen, A. 2008, Journal of Fluid Mechanics, 604, 325 
Gruzinov, A., & Waxman, E. 1999, ApJ, 51 1, 852 
Gruzinov, A. 2008, arXiv:0803.1182 

Kazantsev, A. P. 1968, Soviet Journal of Experimental and Theoretical 

Physics, 26, 1031 
Kulkarni, S. R., et al. 1998, Nature, 395, 663 
Kulsrud, R. M., & Anderson, S. W. 1992, ApJ, 396, 606 
Kulsrud, R. M., Cen, R., Ostriker, J. P, & Ryu, D. 1997, ApJ, 480, 481 



Kulsrud, R. M., & Zweibel, E. G. 2008, Reports on Progress in Physics, 71, 
046901 

Landau, L. D., & Lifshitz, E. M. 1959, Fluid Mechanics (Oxford: Pergamon) 
Medvedev, M. V., & Loeb, A. 1999, ApJ, 526, 697 

Mirabel, I. F, Rodriguez, L. R, Cordier, B., Paul, J., & Lebrun, F 1992, 
Nature, 358, 215 

Nishikawa, K.-L, Hardee, P., Richardson, G., Preece, R., Sol, H., & Fishman, 

G. J. 2003, ApJ, 595, 555 
Pudritz, R. E., & Silk, J. 1989, ApJ, 342, 650 
Rees, M. J., & Setti, G. 1968, Nature, 219, 127 

Schekochihin, A. A., Cowley, S. C, Taylor, S. R, Maron, J. L., & 

McWilliams, J. C. 2004, ApJ, 612, 276 
Soderberg, A. M., et al. 2008, Nature, 453, 469 
Spitkovsky, A. 2008, ApJ, 682, L5 



